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Summary: In this work, the fiber kinking phenomenon, which is known as the failure 
mechanism that takes place when a fiber reinforced polymer is loaded under longitudinal 
compression, is studied. A computational micromechanics model is employed to interrogate 
the assumptions of a recently developed mesoscale continuum damage mechanics (CDM) 
model for fiber kinking based on the deformation gradient decomposition (DGD) and the 
LaRC04 failure criteria. 


1 INTRODUCTION 


There is a need for more accurate mesoscale models for predicting the fiber kinking failure 
mode in fiber reinforced polymer (FRP) laminates for use in progressive damage analysis 
(PDA) codes. One factor limiting the accuracy of predictions by many state-of-the-art PDA 
codes when fiber kinking is active is that most of the physical characteristics of the fiber 
kinking process are ignored. The conventional continuum damage mechanics (CDM) approach 
uses the same phenomenological model used commonly for longitudinal tension for 
longitudinal compression also [1, 2]. The constitutive law typically follows a linear or bilinear 
softening up to the material point is fully damaged and a stress-free state is reached. Such a 
model is not reasonable for compression failure modes in FRP laminates since it is not expected 
that stress will reduce to zero. Instead, it is intuitive that as long as the material remains in 
contact with itself, some residual stress state will exist as damage evolves. 


The fiber kinking theory introduced by Budiansky offers a physically based approach to 
kinkband initiation and propagation [3, 4]. Kinking theory identifies the relevant mechanisms 
in kinkband formation as a combination of an initial fiber misalignment, shear stress-strain 
behavior that is nonlinear, and large fiber rotation. Fiber kinking theory produces the 
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characteristic mechanical response shown in Figure la where, once the compressive strength 
X,, is reached, the stress drops unstably to a plateau crush stress level. The kink band is 
idealized as shown in Figure 1b with a fiber misalignment angle, @, band angle, 6, and band 
width, w; ,. If the shear nonlinearity follows a Ramberg-Osgood behavior, y = (t + at"”)/G, 

Budiansky's theory predicts the strength as, 

G2 
1 i (1) 
G12 Po\ 1 
1/n (412 Po 

Lepyan ( Esl ) 


where, G12, is the in-plane shear modulus, @ and 7 define the nonlinear shear stress-strain 
curve, and @p is the initial fiber misalignment. It is important to notice that this closed-form 
solution assumes f = 0°. 
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Figure 1. The characteristic constitutive response predicted by fiber kinking (a) and the idealization of the kink 
band formation (b) 


Validation of models for fiber kinking is a challenging undertaking due to the instability of 
the process. Most experimental configurations exhibit unstable failure or interaction of kinking 
with other damage mechanisms yielding limited data or challenging cases for model validation. 
In the absence of detailed experimental investigations, computational micromechanics has 
much utility in providing an alternative basis for evaluating assumptions of models derived at 
the mesoscale. Recently, high-fidelity three-dimensional (3-D) computational 
micromechanical models of fiber kinking have been introduced [5, 6] and have shown 
promising as sources of insights into the fiber kinking process. 


In this paper, a computational micromechanical model for fiber kinking is interrogated to 
evaluate the assumptions of a mesoscale model proposed and further developed by Bergan and 
Leone [7]. Some of the details of this constitutive model are discussed in Section 2. In the 
following section, the computational micromechanical model developed by Naya et al. [5] is 
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adapted for use in evaluating the assumptions of the mesoscale model. Finally, in Section 4 a 
direct comparison between the two models is presented. 


2 MESOSCALE CONSTITUTIVE MODEL 


The mesoscale constitutive model represents the fiber kinking phenomenon considering 
geometric and shear nonlinearities. The model includes the kinematics of the fiber kinking 
process by tracking fiber misalignment, ~, throughout loading. The characteristic mechanical 
response shown in Figure la is not prescribed in the model, but is a result of the shear 
nonlinearity and the fiber rotation. No longitudinal compression fracture energy needs to be 
specified. The details of the implementation are presented in detail in [7, 8]. 


The initial misalignment angle, @po, accounts for fiber misalignments and other 
manufacturing defects that contribute to fiber kinking initiation. Rearranging Eq. (1), the initial 
fiber misalignment is, 


a_i 
Qo = (= yr" (2) 
Giz \Xcn ain 


Material models that exhibit strain-softening behavior are susceptible to mesh sensitivity 
when strain localizes. In conventional CDM models, this deficiency is often addressed with 
Bazant’s crack band theory [9] in which the energy dissipated is scaled by the element size. In 
the present model, there is no crack surface on which traction goes to zero and therefore the 
crack band theory is not applicable. Nonetheless, there is an inherent mesh sensitivity since the 
model includes a strain-softening response leading to strain localization in a band of elements 
after the strength is reached. The method used herein is analogous to the strain decomposition 
method by Costa et al. [10] and follows previous work by Bergan and Leone [7]. The model 
has been implemented in Abaqus/Explicit [15] as a user material subroutine (VWUMAT). 


The kink band width, w,y, is assumed to be smaller than the element size such that the 
element can be decomposed into an undamaged material region and a kink band region, as 
shown in Figure 2. In the kink band region, shear nonlinearity is enabled, whereas in the 
undamaged material region the shear response is linear. The deformation gradient 
decomposition (DGD) approach [7, 8] is used to enforce continuity and equilibrium between 
the undamaged and kink band regions. The DGD approach is enabled (and thus the element is 
decomposed) when plastic strain becomes non-negligible. However, kink band width, wzp, 
cannot be predicted by the model and it is required as an input. Values reported in the literature 
range from 50 um to 200 um from experimental observations [11-14] and numerical models 
[5, 6]. 


Miguel Herraez, Andrew C. Bergan, Carlos Gonzalez, Claudio S. Lopes 


Gurrent Reference 


Kink band region Wkb 


Undamaged material region 


X2 Ip 


Figure 2. Schematic representation of the model: decomposition into the kink band and the undamaged region 


3 COMPUTATIONAL MICROMECHANICS MODEL 


In this section, the micromechanical finite element model (FEM) is described including the 
geometry, discretization, and material properties. Subsequently, a calibration study that was 
conducted to demonstrate equivalence between the shear behavior of the mesoscale and 
micromechanical model is described. 


A 3-D single-fiber micromechanical model is used to interrogate assumptions of the 
mesoscale model. The micromechanical FEM is an extension of the single-fiber model 
described in Naya et al. [5] where it was demonstrated that a single-fiber representative volume 
element (RVE) produces strength predictions in good agreement with a multi-fiber RVE. The 
3-D single-fiber model is used here since the model is a good compromise between 
computational expense and accuracy. A trade-study between several modeling strategies 
concluded that a single-fiber and multi-fiber two dimensional models were less accurate than 
3-D single fiber models. Multi-fiber 3-D models were discarded due to convergence difficulties 
and computational expense. 


3.1 Single-fiber 3-D model for fiber kinking 


The model represents a single carbon fiber extruded in the longitudinal z-direction along 
the half wavelength of a sinusoidal curve of length, L, as shown in Figure 3. The initial 
misalignment is geometrically introduced according to, 


y(z) =L ua (1 — cos (x -)) (3) 


1 


y'(2) = gosin(x=) (4) 
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such that the initial misalignment varies along the length of the fiber with the material 
orientation given by Eq. 4, see Figure 3d. 


The fiber diameter is 7.09 um and the fiber volume fraction is 65%. The model is discretized 
using 8-node fully integrated isoparametric elements (C3D8). The in-plane mesh size is around 
1 um, while in the longitudinal direction it is 10 um, and the model length is set to 500 um. 
Periodicity of the mechanical fields is guaranteed by the application of periodic boundary 
conditions (PBC). Gutkin et al. [16] showed that PBC can be applied on single-fiber models 
for longitudinal compressive strength prediction, X¢, at the expense of inducing B = 0°. 
Simulations were conducted in Abaqus/Standard under a dynamic implicit scheme. 


Fibre sections/orientations 


Partitions normal to the axis 
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Figure 3. Illustration of the single-fiber 3-D model (a), detail of the mesh (b), exploded cut view of the model 
(c), and side view with detail of the longitudinal mesh and material orientation (d) 


The material system used in this study is AS4/8552 carbon epoxy. Carbon fibers are 
assumed to behave as linear elastic transversely isotropic solids with the elastic constants 
shown in Table 1. The matrix behavior is represented using the Lubliner damaged/plasticity 
model included in Abaqus [15]. This constitutive equation allows the material to behave as 
quasi-brittle when subjected to dominant tensile stress while it shows elastic-plastic behavior 
under pressure confinement and compressive loads. The tensile response is, therefore, linear 
and elastic with modulus and Poisson ratio, F,, and v,,, until the tensile failure stress, d¢9, is 
reached. Beyond this point, a quasi-brittle softening is induced in the material, with G; being 
the matrix fracture energy. Under uniaxial compression, the response is linear up to the initial 
yield limit, d,9. Then, stress hardening takes place until the ultimate stress value is reached, 
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Ocy. The matrix plasticity/damage model parameters used in the simulations are reported in 
Table 2. Fiber-matrix interface failure is taken into account using a cohesive crack approach. 
To this end, a cohesive interaction between the fiber and matrix surfaces is defined. The 
cohesive interaction is governed by a mixed-mode traction separation law where damage onset 
is controlled by a quadratic stress criterion with normal strength, N, and shear strength, S. 
Additionally, isotropic coulomb friction, €, after cohesive failure is included in the fiber-matrix 
interaction. The interface parameters used in the simulations are provided in Table 3. 


Ey fF Er Vi2F V23f Gir Go3F asf Qf 
(GPa) (GPa) (GPa) (GPa) (107° *G=>) doe ?C-*) 
231.6 12.97 0.3 0.46 11.3 4.45 -0.9 72 


Table 1: Material properties of AS4 carbon fibers [17]. 


Em Vn Am Oto Gt co Ocu 
(GPa) (107° °c~4) (MPa) (J/m?) (MPa) (MPa) 
231.6 0.35 52 121 90 176 180 


Table 2: Parameters for the matrix damage/plasticity model [18]. 


N S Ge Ge E 
(MPa) (MPa) (J/m?) (J/m?) 
57 85 7 81 0.4 


Table 3: Properties of the fiber-matrix interface [18]. 


A previous thermal step was applied to introduce residual thermal stresses appearing due to 
the cooling down process after curing. 


3.2 In-plane shear response 


The mesoscale model utilizes a Ramberg-Osgood shear nonlinearity curve [19]. Ideally, the 
parameters that define the shear response are obtained from a test that isolates the behavior of 
a single ply subjected to large shear deformations. However, in the absence of such test data, 
the ASTM D3518 test of a +45° laminate subjected to tensile loads is used to define the shear 
nonlinearity behavior [20]. The +45° laminate test data smears a wide variety of damage 
mechanisms into a single stress-strain curve, including large fiber rotations and delamination, 
which are not desirable to include in the shear nonlinearity characterization. Nonetheless, given 
the +45° laminate as the source of material input data for the mesoscale model, calibration of 
the micromechanical model was performed so that the model produced an equivalent response 
for an RVE in order to facilitate one-to-one comparison of the two models for the fiber kinking. 


An RVE of a +45° laminate, as shown in Figure 4b, was developed with the parameters and 
modeling approach described above. The dilatancy angle of the matrix, W, and the temperature 
drop, AT, were adjusted to reproduce the experimental shear curves. The final response 
achieved from the numerical model, the experimental curve, and the Ramberg-Osgood curve 
fit are nearly identical as shown in Figure 4a. 
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Figure 4. Shear stress-strain curves including snapshots of the plastic strain in the matrix (a) and RVE of the 
+45° model (b). 


4 COMPARISON OF THE MODELS 


The predictions from the micromechanical model were compared with a single element 
analysis using the mesoscale model to understand the role of simplifying assumptions in the 
latter one. In both analyses, a shortening displacement in the longitudinal direction was 
prescribed. 


The mesoscale finite element is made of one C3D8R element with a uniform edge length of 
0.15 mm. The material properties used are provided in Table 4. 


Ey, Ex2 G42 V2 V23 a n X¢ Yo S, 
(GPa) (GPa) (GPa) (MPa?~" ) (MPa) (MPa) (MPa) 
130.6 9.106 4.82 0.3 0.45 2.86 1071! 6.49 1400 215 71 


Table 4: AS4/8552 material properties for mesoscale model. 


The results in terms of stress-strain curves up to the peak load show excellent agreement 
between the two models for initial misalignments ranging from 0.5 to 4° as shown in Figure 5. 
The CDM strength prediction follows the expression in Eq. 1, according to the nonlinear shear 
response described through the Ramberg-Osgood parameters (G2, 7), @). 


The stress level (crushing stress) maintained after the peak load plotted in Figure 5a is 
similar in both models. The sudden load drop observed in the FEM models following the peak 
load is due to the dynamic nature of the numerical analyses, nevertheless the load stabilizes 
later and shows an asymptotic response regardless of the initial misalignment. The CDM 
analyses encountered some convergence difficulties and so not all cases reached the crush 
stress regime, especially for small initial misalignment angles. 
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Figure 5. Comparison of the micromechanical and mesoscale CDM models under longitudinal compression for 
different initial misalignments: stress-strain (a) and strength-misalignment (b) curves. 


Regarding the kinematics of the fiber kinking process, it is important to differentiate 
between the maximum, @, and the average fiber rotation, @, see Figure 6. The maximum fiber 
rotation takes place in the mid-section of the fiber where @ = max ~(z), while the average or 
global rotation is based upon the deflection of the fiber as, @ = arctan(uy /L). The global 
rotation of the single-fiber model is comparable to the rotation of the CDM model as shown in 
Figure 7a, while the maximum fiber rotation is higher all along the analysis. It is observed that 
sudden fiber rotation occurs right after the peak load is reached, then keeps increasing linearly. 


CDM Single-fiber model 
CDM model 


Figure 6. Interpretation of the fiber rotation kinematics due to fiber kinking 


From previous experimental studies [13, 21—23], it is reported that once the kink band is 
generated, it propagates through the specimen at some inclination (kink angle, 10° < B < 30°) 
with an approximately constant width. However, in the single-fiber model kink band 
broadening is observed due to the periodic constraints imposed in the model, as shown in 
Figure 7b. For this reason, kink band width increases progressively and the representative value 
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for Wxp is the one observed right after load drop, which is around 100 wm ~ 15 d;. The kink 
band width was one of the input parameters of the mesomechanical model and w,;, = 100 um 
was selected for this analysis. 
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Figure 7. Comparison of the kinematics in the two models: fiber rotation (a) and kink band width (b). Markers 
represent the point when the peak load is reached in the FEM model: triangles stand for @ and squares for @. 


The effect of fiber-matrix friction coefficient was analyzed by conducting the same 
simulations with a frictionless interface (§ = 0). The results of this analyses are shown in 
Figure 8. Although the compressive strength was barely affected, it was observed that there 
was a drop of approximately 30% in the crushing stress. This drop in stress is explained by the 
more localized kink band (w;, = 70 um — 80 um), whose local rotation, @, increases up to 
30% compared to the reference case with friction. Budiansky et al. [4] related the residual 
crushing stress, 0,,, with the in-plane shear strength of the material, S;, and the kink band 


angle, f, as: 


2s 
Per = Sin 2B 


(5) 


The single-fiber model does not represent £6 due to the periodic boundary conditions. 
Nevertheless, experimental observations from Vogler et al. [24] showed that f is proportional 
to the fibers rotation, g. Thus, larger fiber rotation induces higher kink band angle sustaining 


lower crushing stress levels. 
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Figure 8. Results comparing the effect of the fiber-matrix friction on the fiber kinking mechanism: stress-strain 
curves (a) and rotation-strain curves (b) for different initial misalignments. Markers represent the point when the 
peak load is reached in the FEM model: triangles stand for @ and squares for @. 


5 CONCLUSIONS 


In this study a computational micromechanics model is used to assess the assumptions used 
in the derivation of a continuum damage mechanics model considering fiber kinking failure 
mechanism. The continuum damage mechanics (CDM) model employs the deformation 
gradient decomposition (DGD) strategy and is able to capture the main features of fiber kinking 
such as sudden load drop, kink band formation and rotation, and holds a residual stress level 
(crushing). 


Very good agreement was found between both models not only in terms of compressive 
strength prediction, but also regarding the kinematics of fiber kinking (rotation and band width) 
and crushing stress. 
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